{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DoWhy example on Twins dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we study the twins dataset as studied by <a href=\"https://arxiv.org/pdf/1705.08821.pdf\" target=\"_blank\">Louizos et al</a>. We focus on twins which are the same sex and weigh less than 2kgs. The treatment t = 1 is being born the heavier twin and the outcome is mortality of each of the twins in their first year of life.The confounding variable taken is 'gestat10', the number of gestational weeks prior to birth, as it is highly correlated with the outcome. The results using the methods below are in coherence with those obtained in the paper. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import dowhy\n",
    "from dowhy import CausalModel\n",
    "from dowhy import causal_estimators\n",
    "\n",
    "# Config dict to set the logging level\n",
    "import logging.config\n",
    "DEFAULT_LOGGING = {\n",
    "    'version': 1,\n",
    "    'disable_existing_loggers': False,\n",
    "    'loggers': {\n",
    "        '': {\n",
    "            'level': 'WARN',\n",
    "        },\n",
    "    }\n",
    "}\n",
    "logging.config.dictConfig(DEFAULT_LOGGING)\n",
    "# Disabling warnings output\n",
    "import warnings\n",
    "from sklearn.exceptions import DataConversionWarning\n",
    "warnings.filterwarnings(action='ignore', category=DataConversionWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font size=\"4\">Load the Data</font>**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data loading process involves combining the covariates, treatment and outcome, and resolving the pair property in the data. Since there are entries for both the twins, their mortalities can be treated as two potential outcomes. The treatment is given in terms of weights of the twins.Therefore, to get a binary treatment, each child's information is added in a separate row instead of both's information being condensed in a single row as in the original data source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The covariates data has 46 features\n",
    "x = pd.read_csv(\"https://raw.githubusercontent.com/AMLab-Amsterdam/CEVAE/master/datasets/TWINS/twin_pairs_X_3years_samesex.csv\")\n",
    "\n",
    "#The outcome data contains mortality of the lighter and heavier twin\n",
    "y = pd.read_csv(\"https://raw.githubusercontent.com/AMLab-Amsterdam/CEVAE/master/datasets/TWINS/twin_pairs_Y_3years_samesex.csv\")\n",
    "\n",
    "#The treatment data contains weight in grams of both the twins\n",
    "t = pd.read_csv(\"https://raw.githubusercontent.com/AMLab-Amsterdam/CEVAE/master/datasets/TWINS/twin_pairs_T_3years_samesex.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#_0 denotes features specific to the lighter twin and _1 denotes features specific to the heavier twin\n",
    "lighter_columns = ['pldel', 'birattnd', 'brstate', 'stoccfipb', 'mager8',\n",
    "       'ormoth', 'mrace', 'meduc6', 'dmar', 'mplbir', 'mpre5', 'adequacy',\n",
    "       'orfath', 'frace', 'birmon', 'gestat10', 'csex', 'anemia', 'cardiac',\n",
    "       'lung', 'diabetes', 'herpes', 'hydra', 'hemo', 'chyper', 'phyper',\n",
    "       'eclamp', 'incervix', 'pre4000', 'preterm', 'renal', 'rh', 'uterine',\n",
    "       'othermr', 'tobacco', 'alcohol', 'cigar6', 'drink5', 'crace',\n",
    "       'data_year', 'nprevistq', 'dfageq', 'feduc6', 'infant_id_0',\n",
    "       'dlivord_min', 'dtotord_min', 'bord_0',\n",
    "       'brstate_reg', 'stoccfipb_reg', 'mplbir_reg']\n",
    "heavier_columns = [ 'pldel', 'birattnd', 'brstate', 'stoccfipb', 'mager8',\n",
    "       'ormoth', 'mrace', 'meduc6', 'dmar', 'mplbir', 'mpre5', 'adequacy',\n",
    "       'orfath', 'frace', 'birmon', 'gestat10', 'csex', 'anemia', 'cardiac',\n",
    "       'lung', 'diabetes', 'herpes', 'hydra', 'hemo', 'chyper', 'phyper',\n",
    "       'eclamp', 'incervix', 'pre4000', 'preterm', 'renal', 'rh', 'uterine',\n",
    "       'othermr', 'tobacco', 'alcohol', 'cigar6', 'drink5', 'crace',\n",
    "       'data_year', 'nprevistq', 'dfageq', 'feduc6',\n",
    "       'infant_id_1', 'dlivord_min', 'dtotord_min', 'bord_1',\n",
    "       'brstate_reg', 'stoccfipb_reg', 'mplbir_reg']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Since data has pair property,processing the data to get separate row for each twin so that each child can be treated as an instance\n",
    "data = []\n",
    "\n",
    "for i in range(len(t.values)):\n",
    "    \n",
    "    #select only if both <=2kg\n",
    "    if t.iloc[i].values[1]>=2000 or t.iloc[i].values[2]>=2000:\n",
    "        continue\n",
    "    \n",
    "    this_instance_lighter = list(x.iloc[i][lighter_columns].values)\n",
    "    this_instance_heavier = list(x.iloc[i][heavier_columns].values)\n",
    "    \n",
    "    #adding weight\n",
    "    this_instance_lighter.append(t.iloc[i].values[1])\n",
    "    this_instance_heavier.append(t.iloc[i].values[2])\n",
    "    \n",
    "    #adding treatment, is_heavier\n",
    "    this_instance_lighter.append(0)\n",
    "    this_instance_heavier.append(1)\n",
    "    \n",
    "    #adding the outcome\n",
    "    this_instance_lighter.append(y.iloc[i].values[1])\n",
    "    this_instance_heavier.append(y.iloc[i].values[2])\n",
    "    data.append(this_instance_lighter)\n",
    "    data.append(this_instance_heavier)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = [ 'pldel', 'birattnd', 'brstate', 'stoccfipb', 'mager8',\n",
    "       'ormoth', 'mrace', 'meduc6', 'dmar', 'mplbir', 'mpre5', 'adequacy',\n",
    "       'orfath', 'frace', 'birmon', 'gestat10', 'csex', 'anemia', 'cardiac',\n",
    "       'lung', 'diabetes', 'herpes', 'hydra', 'hemo', 'chyper', 'phyper',\n",
    "       'eclamp', 'incervix', 'pre4000', 'preterm', 'renal', 'rh', 'uterine',\n",
    "       'othermr', 'tobacco', 'alcohol', 'cigar6', 'drink5', 'crace',\n",
    "       'data_year', 'nprevistq', 'dfageq', 'feduc6',\n",
    "       'infant_id', 'dlivord_min', 'dtotord_min', 'bord',\n",
    "       'brstate_reg', 'stoccfipb_reg', 'mplbir_reg','wt','treatment','outcome']\n",
    "df = pd.DataFrame(columns=cols,data=data)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.astype({\"treatment\":'bool'}, copy=False) #explicitly assigning treatment column as boolean \n",
    "\n",
    "df.fillna(value=df.mean(),inplace=True)    #filling the missing values\n",
    "df.fillna(value=df.mode().loc[0],inplace=True)\n",
    "\n",
    "data_1 = df[df[\"treatment\"]==1]\n",
    "data_0 = df[df[\"treatment\"]==0]\n",
    "print(np.mean(data_1[\"outcome\"]))\n",
    "print(np.mean(data_0[\"outcome\"]))\n",
    "print(\"ATE\", np.mean(data_1[\"outcome\"])- np.mean(data_0[\"outcome\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font size=\"4\">1. Model</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#The causal model has \"treatment = is_heavier\", \"outcome = mortality\" and \"gestat10 = gestational weeks before birth\"\n",
    "model=CausalModel(\n",
    "        data = df,\n",
    "        treatment='treatment',\n",
    "        outcome='outcome',\n",
    "        common_causes='gestat10'\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font size=\"4\">2. Identify</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)\n",
    "print(identified_estimand)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font size=\"4\">3. Estimate Using Various Methods</font>**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font size=\"3\">3.1 Using Linear Regression</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "estimate = model.estimate_effect(identified_estimand,\n",
    "        method_name=\"backdoor.linear_regression\", test_significance=True\n",
    ")\n",
    "\n",
    "print(estimate)\n",
    "print(\"ATE\", np.mean(data_1[\"outcome\"])- np.mean(data_0[\"outcome\"]))\n",
    "print(\"Causal Estimate is \" + str(estimate.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font size=\"3\">3.2 Using Propensity Score Matching</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "estimate = model.estimate_effect(identified_estimand,\n",
    "        method_name=\"backdoor.propensity_score_matching\"\n",
    ")\n",
    "\n",
    "print(\"Causal Estimate is \" + str(estimate.value))\n",
    "\n",
    "print(\"ATE\", np.mean(data_1[\"outcome\"])- np.mean(data_0[\"outcome\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font size=\"4\">4. Refute</font>**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font size=\"3\">4.1 Adding a random cause</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "refute_results=model.refute_estimate(identified_estimand, estimate,\n",
    "        method_name=\"random_common_cause\")\n",
    "print(refute_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font size=\"3\">4.2 Using a placebo treatment</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_placebo=model.refute_estimate(identified_estimand, estimate,\n",
    "        method_name=\"placebo_treatment_refuter\", placebo_type=\"permute\",\n",
    "        num_simulations=20) \n",
    "print(res_placebo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<font size=\"3\">4.3 Using a data subset refuter</font>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_subset=model.refute_estimate(identified_estimand, estimate,\n",
    "        method_name=\"data_subset_refuter\", subset_fraction=0.9,\n",
    "        num_simulations=20)\n",
    "print(res_subset)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
